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ABSTRACT 

We present a general method for determining the masses and orbital parameters 
of binary millisecond pulsars with long orbital periods (-P rb 3> ly r )) using timing 
data in the form of pulse frequency derivatives. Our method can be used even when 
the available timing data cover only a small fraction of an orbit, but it requires 
high-precision measurements of up to five successive derivatives of the pulse frequency. 
With five derivatives a complete determination of the mass and orbital parameters 
is in principle possible (up to the usual inclination factor sini). With less than five 
derivatives only partial information can be obtained, but significant constraints can 
sometimes be placed on, e.g., the mass of the companion. We apply our method to 
analyze the properties of the second companion in the PSR B1620-26 triple system. 
We use the latest timing data for this system, including a recent detection of the fourth 
time derivative of the pulse frequency, to constrain the mass and orbital parameters of 
the second companion. We find that all possible solutions have a mass mi in the range 
2.4 x 10 -4 Mq < ni2sini2 < 1.2 x 10~ 2 M Q , i.e., almost certainly excluding a second 
companion of stellar mass and suggesting instead that the system contains a planet or a 
brown dwarf. To further constrain this system we have used preliminary measurements 
of the secular perturbations of the inner binary. Using Monte-Carlo realizations of the 
triple configuration in three dimensions we find the most probable value of mi to be 
0.01 ± 0.005 Mq, corresponding to a distance of 38 ± 6AU from the center of mass 
of the inner binary (the errors indicate 80% confidence intervals). We also apply our 
method to analyze the planetary system around PSR B1257+12, where a distant, giant 
planet may be present in addition to the three well-established Earth-mass planets. We 
find that the simplest interpretation of the frequency derivatives implies the presence 
of a fourth planet with a mass of ~ 100 M ffi in a circular orbit of radius ~ 40 AU. 
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1. Introduction 

The traditional method for obtaining masses and orbital parameters of binary pulsars consists 
of fitting Keplerian or Post-Newtonian models to timing data covering at least one complete 
orbit. For wide-orbit binary pulsars, with orbital periods longer than a few years to decades, 
fitting a complete orbit may not be possible. For such systems, we present a method for obtaining 
the masses and orbital parameters using measured values of the successive time derivatives of 
the pulse frequency (/, /, /, etc.). Given the high precision of millisecond pulsar timing it is 
sometimes possible to measure these frequency derivatives up to high order with only a few years 
of observations. We show below that with five derivatives a complete solution may be obtained, 
with the orbital parameters and companion mass fully determined up to the usual unknown 
inclination factor sini. This solution is constructed under the assumption that all frequency 
derivatives are dynamically induced rather than being intrinsic to the pulsar spin-down. This is 
a crucial assumption which may not always be justified. Its validity, and the effects of relaxing 
it, must be examined for each particular application. If only three or four derivatives have been 
measured, significant constraints can still be placed on the parameters of the system. With just 
three dynamically-induced derivatives and the assumption of a circular orbit (often justified, e.g., 
for a planet) a complete solution can again be obtained. Our method is not a substitute for 
the standard fitting procedure to data covering a complete orbit. Instead, it provides a way of 
obtaining the orbital parameters of wide-orbit binaries which cannot be observed over a complete 
orbit due to their very long orbital periods. The method can only be successfully applied to 
binaries containing fast millisecond pulsars with low timing noise, in which one can reasonably 
expect the dynamically-induced pulse frequency derivatives to dominate over intrinsic changes. 

Applications of our method to two systems, PSR B1620-26 and PSR B1257+12, will be 
presented in this paper. Both of these systems are thought to contain planetary-mass companions 
to the millisecond pulsar (Wolszczan 1994; Arzoumanian et al. 1996). Their properties are 
of great interest for our understanding of planet formation outside the solar system. This is 
particularly true in light of the many recent detections of extrasolar planets around nearby 
stars, which show a great diversity of properties (Mayor h Queloz 1995; Butler & Marcy 1996; 
Marcy & Butler 1996). PSR B1620— 26, in the globular cluster M4, is in a hierarchical triple 
configuration. Previously available timing data allowed a second companion mass anywhere in the 
range ~ 10 _3 -1 M (Michel 1994; Sigurdsson 1995), i.e., including the possibility of a Jupiter- type 
planet. PSR B1257+12 has three confirmed low-mass planets, and there is recent evidence for a 
fourth more massive one (Wolszczan 1996). 

Our paper is organized as follows. In §2, we describe our general method for obtaining the 
companion mass and the orbital parameters using measured pulse frequency derivatives. In §3, 
we apply our method to PSR B1620— 26. We also incorporate into our analysis preliminary 
measurements of secular perturbations of the inner binary by performing Monte-Carlo simulations 
with the undetermined parameters in the system to obtain the most probable mass for the second 
companion. In §4, we apply our method to check whether the observed frequency derivatives of 
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PSR 1257+12 are consistent with the possible existence of a fourth object in the planetary system. 

2. Inverting Frequency Derivative Data 

In the standard method for determining the parameters of binary pulsars, the data (consisting 
of radio pulse arrival times) are fitted to the predictions of a Keplerian or Post-Newtonian model 
of the binary orbit. To obtain a reliable fit, one usually needs timing data covering at least a 
few complete orbital periods. However, even if the data do not cover a full period, it is still 
possible to determine, at least approximately, the companion mass and orbital parameters of the 
system if sufficiently accurate timing data are available. The method we develop here uses time 
derivatives of the pulse frequency (basically the coefficients in a Taylor expansion of the pulse 
frequency around a particular epoch). Pulse frequency derivatives are a convenient way in which 
radio astronomers can present the results of their observations when a clear periodicity cannot be 
recognized in the timing data. 

2.1. General Formulation 

Assuming that the pulsar mass (mi ~ 1.4M ) is known, there are five parameters that can 
in principle be determined using our method: the mass of the companion m-2 (up to the unknown 
inclination angle 12), the semi-major axis 02, the eccentricity e2, the longitude of pericenter u>2 
(measured from the ascending node) and the longitude A2 of the companion at the reference 
epoch (measured from pericenter). The inclination angle 12 is the angle between the normal to 
the orbital plane and the line of sight; this angle cannot be determined directly from the timing 
data. Here and throughout this paper a subscript 1 refers to the pulsar, a subscript 2 refers to the 
companion, and all orbital elements correspond to the motion around the center of mass of the 
system (e.g., the distance between the pulsar and its companion is r\ + r2). 

If frequency derivatives up to the fifth order are measured, all five parameters (r«2, 02, e2, 
UJ2, A2) can in principle be determined. For small values of 1712/ mi the usual combination 1712 sin 12 
can be obtained (see the discussion for PSR B1620— 26 in §3). For larger companion masses, the 
dependence of the solutions on 12 is more complicated and one needs to adopt a particular value 
of sin 22 (in practice sin 22 ^ 1 for a random orientation) in order to solve explicitly for 7712- 

Our method assumes that the measured frequency derivatives are purely those induced by 
the motion of the pulsar (acceleration, jerk, and higher derivatives) around the binary center of 
mass. This requires correcting the measurements for other possible kinematic effects (see, e.g., 
Camilo et al. 1994). More importantly, one must assume that any intrinsic contribution from the 
pulsar spin-down can be neglected, or determined to some extent from the known properties of 
other millisecond pulsars and subtracted from the measured values. In particular, the observed 
first frequency derivative / is in general determined by a combination of acceleration and intrinsic 
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spin-down of the pulsar, 

fobs = f'mt ~\~ face- (1) 

It is not possible to measure the two components separately in general. For some systems, 
however, it is reasonable to assume that |/i nt | <C |/ aC c| if |/obs| is large. For example, the observed 
/obs may be positive, a clear sign that it is determined predominantly by acceleration (as in the 
M15 pulsars PSR 2127+11 A and D; Wolszczan et al. 1989). The expected value of /j nt may also 
be estimated from the pulse frequency / and from the assumption that the timing age r = f/(2f) 
of a millisecond pulsar should satisfy r ^ 10 9 yr. Indeed most millisecond pulsars with reliably 
measured timing ages appear to satisfy this property (e.g., Phinney & Kulkarni 1994). Note that 
the true ages of some millisecond pulsars may be considerably smaller (cf. Lorimer et al. 1995), 
but this does not affect our argument. Similarly, the expected value of /i n t can be estimated 
from the predicted level of timing noise for the pulsar, which, although very small for millisecond 
pulsars, may also affect the interpretation of / Q b s (see, e.g., Arzoumanian et al. 1994; Kaspi et al. 
1994). Intrinsic higher derivatives (/, etc.) are normally not measurable for millisecond pulsars, 
and therefore we can always assume safely that any measured values are dynamically induced. 

After subtraction of known kinematic and intrinsic contributions, we can write the time 
derivatives of the pulse frequency at a particular reference epoch as 

/ = -/^, (2) 

/ = -/— + £ etc., 
c / 

where c is the speed of light, a is the acceleration of the pulsar, and n is a unit vector in the 
direction of the line of sight, and a dot indicates a time derivative. The second term in the 
expression for / is smaller than the first by a factor ~ v/c, where v is the orbital velocity of the 
neutron star, i.e., f 2 /f <C |/|. For /, we have \ff/f \ <C |/|, etc., so that all similar terms can be 
neglected in taking higher and higher derivatives. Therefore, we can write the first five frequency 
derivatives simply as 

/ = 

/ = /"• (3) 



7 = -/— . 

c 

Equations (3) form a system of five nonlinear algebraic equations with five unknowns (rri2, 02, &ii 
A2, and CJ2) which must be solved numerically. It is straightforward but tedious to write down 
explicitly the right-hand sides of these equations in terms of the five unknowns, and we will omit 
their explicit forms in the general case. Some of the steps involved, however, are described below 
in §2.2. 
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2.2. Solution with Four Derivatives 

When four derivatives (/ through /) are known, one can obtain a one-parameter family of 
solutions. In practice, the problem can be reduced to solving a system of three nonlinear equations 
as follows. If mi is the mass of the pulsar, and ?rt2 is the mass of the companion, the acceleration 
of the pulsar in its motion around the center of mass of the binary is of magnitude a = k/rf, 
where r\ is the distance from the pulsar to the center of mass and 

k = G ( m } (4) 

For an elliptic Keplerian orbit of semimajor axis a\ and eccentricity e\ = e 2 , the distance r\ is 
given by 

11 A 

— = -(1 + e 2 cosAi) = — , (5) 
r\ h h 

where h = ai(l — e|). Equations (3) for the pulse frequency derivatives can then be written 



,a 



f = — f- sin(Ai + ll>\) sini2 = —f'KA sin(Ai + cji), (6) 

/ = -fKBM, (7) 

/ = -fKCXl (8) 

and '/ = -fKDXl (9) 



where we have defined 



A = l + e2CosAi, 

B = 2AA'sm.{\i+uJi)+A 2 cos{\i+uJi) 1 

C = Bf+^f, (10) 

D = C' + 4CA ' 



and K 



A 
k sin i 2 



h?c ' 

and a prime indicates a derivative with respect to Ai. 

Using equation (6), we can rewrite equations (7)-(9) as 



f = EM (id 

f = 9M (12) 

1 yl 2 sin(Ai+ Wl )' 1 ' 

and J = A9 D jf lf -. (13) 

J ^ 2 sin(Ai +wi) v ; 
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This is a nonlinear system of three equations with four unknowns- e 2 , Ai = A2, 0J\ = ^2 + 180°, 
and Ai. Assuming a value for one of them we can solve for the remaining parameters. 

We have chosen to use the eccentricity e2 as our free parameter. For an assumed value of e2, 
we solve equations (11)— (13) for Ai, Ai, and u\ using the Newton-Raphson method (see e.g., Press 
et al. 1992). This method requires an initial guess for the unknown parameters, which is then 
successively improved until convergence to an actual solution is obtained. One must be careful to 
experiment with many different initial guesses, since nonlinear systems like this one often have 
multiple branches of solutions. In addition, symmetries must also be taken into account. Here 
physically equivalent solutions are obtained if the direction of motion is reversed, and the signs of 
u)\ and Ai are also reversed. This gives a different but equivalent orientation of the system. Once 
Ai, Ai, and u\ are known, conservation of angular momentum together with equation (6) gives 

A = Ai (14) 

_ k f c , , 

h 2 fA 2 sm(X 1 +w 1 )sinz 2 



Using equations (14) and (15), we can calculate k and h, 

fcA 2 
/sin 72 sin AiA^ 



fcA 2 

h = - — , . (16) 



and k = -(—-A—) 3 (*L). (17) 

\fsmi 2 sm\ 1 J \\fj y J 

It is now straightforward to obtain 777,2 (from k, assuming that the pulsar mass mi is known) and 
the semi-major axis a>2 = (mi/m2)/i/(l — e|). When m2 mi, we can directly obtain 7772 sin 12 as 
follows. For small 7772, k pa Gmf/mf. Therefore 7772 ~ (km\/G) 1 / 3 and equation (17) gives 

fc fm 2 A 2 \ 1/3 , , 

m 2 sim2^—. — — — ^r— . (18) 
/sinAi \GXf J 

For larger values of 7772/7771, one needs to assume a value for sin 72 and solve equation (4) 
numerically for 7772. 



2.3. Solution with Three Derivatives 

When only three derivatives are known, one can assume values for two parameters and solve 
for the remaining three parameters in the same way as above. Alternatively, in the special case of 
a circular orbit (e 2 = 0), the system can be solved completely (to within sin 72) using only three 
derivatives. In that case, we have 

• -fksini 2 . . , , . 

/ = 2 sm (Ai), (19) 

CLi c 
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and 



/ 

7 



-fksmi 2 , 

2 cos(Ai)Ai, 

a\ c 

fksmi2 



sin(Ai)Af . 



at c 



(20) 
(21) 



Note that here Ai is the longitude of the pulsar measured from the ascending node (not from 
pericenter as before). 



Using equation (19) to eliminate — /fesin^/af c, we get 



/ 



and 



^ sin(Ai) 
f = /At. 



cos(Ai)Ai, 



This gives 



A? = -(///), 



and 



Ai = arctan 



f -f 



1/2- 



./ V / / 

Using equation (19) and conservation of angular momentum, we have 



(22) 
(23) 

(24) 
(25) 



A: 



-fc 



f sini2 sin(Ai^ 



and 



A: 



- i 2 
3 _ A i- 



(26) 
(27) 



Dividing equation (26) by equation (27), and then using equations (24), (25) and (27), we get 

f 2 c 



and 



ai 



A; 



/ / sini 2 sin(Ai) 
\ / sin %2 sin(Ai 



(28) 
(29) 



where Ai is given by equation (25). For -C mi, A; w Gm%/m\ . Then, for given values of the 
frequency derivatives and the pulsar mass mi, m 2 can be calculated explicitly using equation (29). 
We find 



7722 sin 22 



' 2 \ V3 



fc 



2/3 



(30) 



where Ai is given by equation (25). Finally we can calculate a2 = (m 1 /m 2 )ai, where a\ is given 
by equation (28). 



-8- 



2.4. Solution with Two Derivatives 

When only two derivatives are available, one can obtain a one-parameter family of solutions, 
assuming again that e<i ss 0. We can use / (or equivalently Ai, using eq. [24]) as the free parameter 
and then use the equations of §2.3 to construct a one-parameter family of solutions explicitly. 

3. Application to the PSR B1620-26 Triple System 

The millisecond pulsar PSR B1620— 26, in the globular cluster M4, has a low-mass binary 
companion (probably a white dwarf of mass m c ss 0.3 M Q for a pulsar mass m p = 1.35 M ) in a 
191 day low-eccentricity orbit (Lyne et al. 1988; McKenna & Lyne 1988). The unusually large 
frequency second and third derivatives indicate the presence of a second companion around the 
inner binary, forming a hierarchical triple configuration (Backer 1992; Backer, Foster, & Sallmen 
1993; Thorsett, Arzoumanian, & Taylor 1993). Such hierarchical triple systems are expected 
to be produced quite easily in dense globular clusters through dynamical interactions between 
binaries. In a typical interaction, one star would be ejected, leaving the other three in a stable 
triple system (Rasio, McMillan, & Hut 1995; Sigurdsson 1995). Previous calculations (done using 
frequency derivatives up to the third order) showed that the mass of the second companion 
could be anywhere from ~ 10" 3 M to ~ 1M Q (Michel 1994; Sigurdsson 1995). Recently, a 
measurement has been made of the fourth derivative of the pulse frequency, along with preliminary 
measurements of secular changes of the inner binary parameters due to the perturbation of the 
second companion (Arzoumanian et al. 1996). These include a precession of the inner binary 
orbital plane (measured as a change in the projected semi- major axis of the binary), and possible 
changes in the eccentricity and longitude of periastron. 

3.1. Modeling the Frequency Derivatives 

In this section, we apply the method described in §2.2 to analyze the properties of the second 
companion in the PSR B1620— 26 triple system. Since the orbital period of the second companion 
is much longer than that of the inner binary (for all solutions obtained below), we treat the inner 
binary as a single object. Keeping the same notation as before, we let mi = m p + m c be the mass 
of the inner binary pulsar, with m p the mass of the neutron star and m c the mass of the (inner) 
companion, and we denote by mi the mass of the second companion (to be determined). As in §2 
the orbital parameters (A2, ui2, &2, «2 and 12) refer to the orbit of the (second) companion with 
respect to the center of mass of the system (here the entire triple). However, a subscript 1 for the 
orbital elements refers to the orbit of the inner binary. The results presented in this section are 
all for mi = 1.7 Mq (assuming m p = 1.4 M Q , m c = 0.3 M ; cf. Thorsett et al. 1993). However, 
we have have checked that they are not very sensitive to small changes in the value of m\. In 
particular the companion mass mi varies only as ~ m/ (cf. eq. [4]). 
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We use the latest available values of the pulse frequency derivatives (Arzoumanian & Thorsett 
1996) for the epoch MJD 48725.0: 



/ 


= 90.2873320054(1) s" 1 


/ 


= -5.4702(7) x 10~ 15 s 


/ 


= 1.929(8) x 10~ 23 s- 3 


7 


= 8(1) x 10~ 33 s" 4 


J 


= -2.1(6) x KT 40 s- 5 



These values take into account a (very precise) Keplerian model of the inner orbit. The corrections 
to / due to proper motion are negligible for this pulsar. The frequency derivatives should 
therefore reflect the residual motion of the pulsar caused by the presence (unmodeled) of a second 
companion. However, as discussed in §2.1, the observed first derivative / Q bs can be a combination 
of the intrinsic spin-down of the pulsar and the acceleration due to the second companion. Since 
/obs is negative, and /i nt is always negative, / acc can in principle be either positive or negative. If 
it is negative, then its magnitude must be <| / b s |- If it is positive, its magnitude can in principle 
be larger. However, it cannot be much larger than | / Q b s | since this would imply a very large 
intrinsic spin-down rate, and a short characteristic age r = — //(2/), which is not expected for a 
millisecond pulsar (cf. Thorsett et al. 1993 and §2.1). In practice, we find that varying /acc/ /obs in 
the entire range —1.0 to +1.0 does not affect our solutions significantly (see below). The expected 
value of / due to timing noise can be estimated using, e.g., Figure 1 of Arzoumanian et al. (1994), 
which gives the timing noise parameter As = log(|/|10 24 /6/) as a function of period derivative. 
This gives an upper limit on the contribution to / due to intrinsic timing noise of 3 x 10~ 24 s -3 , 
which is an order of magnitude smaller than |/ bs|- The same conclusion is reached if we consider 
for comparison PSR B1855+09, which has a comparable spin rate (/ = 186 s -1 ) but a frequency 
second derivative | jobs I < 2 x 10 _27 s~ 3 (3a) i.e., at least four orders of magnitude smaller than 
/obs in PSR B 1620-26 (Kaspi et al. 1994). 

Figure 1 illustrates our "standard solution" , obtained under the assumption that / acc = / D b s 
and using the current best-fit value for the fourth derivative, f m = —2.1 x 10 -40 s -5 . Following the 
method of §2.2, we use as the free parameter. We see that there are no solutions for e2 5s 0.1. 
Hence a nearly circular orbit is ruled out. For 0.1 ^ ei ^ 0.3 there are two solutions for each value 
of the eccentricity, and hence two possible values of ni2. 

In the first branch of solutions, mi approaches zero as the eccentricity approaches the 
value e2 ~ 0.33. However for very small mi the second companion gets close enough to the 
inner binary to make the triple configuration dynamically unstable. The stability of the triple 
system can be checked using an approximate criterion of the form Y > Y m [ n (for stability), where 
Y = [(1 — ei)(a\ + ai)]/[(a p + a c )(l + ei)] is the ratio of the outer pericenter separation (from 
the center of mass of the inner binary to the second companion) to the apastron separation of the 
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inner binary. Here we have used the results of Eggleton & Kiseleva (1995), who give 

3.7 2.2 1.4 g^J - 1 
Ymia ~ 1 + 7?3 +7— T73+773 773-7- (32) 

Vout 1 T Vout Sin Vout T 1 

Here c/i n = m p /m c and (fout = 7711/777,2 are the inner and outer mass ratios. Using the values 
m p ~ 1.4M0, m c ~ O.3M0, and ni2 ~ 1O _5 M0, we get Y m \ n ~ 2. Thus we require that Y > 2 
for a dynamically stable solution. This gives us a lower limit on the second companion mass 
771-2 ^ 3 x 10~ 5 M w 10 M®. In addition, we can rule out solutions with orbital periods P2 & 14 yr 
since this would have been detected already in the timing residuals (S. E. Thorsett, personal 
communication). This gives us a somewhat stricter lower limit of 7772 ^ 2.4 x 10~ 4 M « 80M e . 

For the second branch of solutions, 7772 increases monotonically from ~ 10~ 3 M to 
1.2 x 10~ 2 M (i.e., Jupiter to brown-dwarf masses) as e2 increases from ~ 0.1 to 1. As e2 
approaches 1.0, the mass 7772 remains bound. So even though the solutions with e2 — > 1 are a 
priori unlikely, they provide a strict upper limit 771-2 < 1-2 x 1O~ 2 M0. Since all our solutions have 
777-2 *C 771-1, our method gives in fact the product 7772 sin 22 (cf. §2.2), which we show in Figure 1. 

The effect of varying / acc //obs in the range 0-1.0 is shown in Figure 2. We see that 
smaller values of / acc //obs give slightly smaller mass solutions, but the overall mass range is not 
significantly affected. In particular, we see that stellar-mass solutions for 777,2 are not allowed even 
if the observed value of / is assumed to be mostly intrinsic. Negative values of / aC c//obs give 
similar results. 

We see that our standard solution excludes the stellar mass range (7772 ^ 0.1 M ) for 
the second companion, and this is rather surprising. A stellar mass would provide a natural 
explanation for the anomalously high eccentricity of the inner binary (ei ~ 0.03) in terms of 
secular perturbations (Rasio 1994). In addition, a stellar mass would also be consistent with 
a preliminary identification of an optical counterpart for the system (Bailyn et al. 1994). We 
have seen already that assuming a different value for / acc does not change our conclusions. 
Alternatively, we can try to vary the value of / within its fairly large error bar. However, we find 
that varying / within its formal la error bar still does not produce significant changes in 7112. In 
order to obtain stellar-mass solutions, we find that we must vary / within a larger 4<r error bar 
around the best fit value f m given above. Note that this allows for a change of sign in the actual 
value of /. The results are illustrated in Figure 3. Here we assume that sin 72 = 1- We find that 
stellar-mass solutions (with 7712 ^ 0.1M Q ) are possible only if —0.1 <;/ / f m ^ 0.1, i.e., / must 
be more than 3a away from its present best-fit value (assuming sin 7-2 ~ 1). 



3.2. Secular Perturbations of the Inner Binary 



We can further constrain the system by considering secular perturbations in the orbital 
elements of the inner binary. Preliminary measurements have been made of the perturbations 
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in oji, ei, and xi (Arzoumanian al. 1996; Arzoumanian & Thorsett 1996). We use these 
measurements to further constrain the system by requiring that all our solutions be consistent 
with these secular perturbations. 

Since the period P2 of the second companion is much larger than the period of the inner 
binary, we can calculate the secular perturbations assuming that the second companion has a fixed 
position in space with respect to the inner binary. Let r\2, 02 and 4>2 be the fixed spherical polar 
coordinates of the second companion, with the origin at the center of mass of the inner binary, 4>2 
measured from pericenter in the orbital plane, and 62 such that sin 62 = 1 for the coplanar case. 
Then the averaged perturbation rates are given by (Rasio 1994) 

= ^[ sin 2 (0 2 )(5cos 2 (0 2 )-l)-l] (33) 
— 15tt 

ei = — — r/eisin 2 (6>2)sin(2</> 2 ) (34) 
3tt 

h = 7-7^77 sin(26»2)cos(cJi + 4> 2 ) (35) 
lf\ 

where Pi = 16540653s is the period of the inner binary, 77 = (m 2 /mi)[(a p + a c )/Vi 2 ] 3 , a c is the 
semi-major axis of the inner binary companion, and a p is the semi-major axis of the pulsar (with 
respect to the center of mass of the inner binary). The projected semi-major axis of the pulsar is 
x p = (l/c)a p sinii, and therefore 

x p = -a p cos (36) 
Note that there is no secular perturbation of the semi-major axis (a p = 0). 

The present measured values of the perturbations are (Arzoumanian & Thorsett 1996) 

wi = (-2.0 ±2.1) x 10- 4o yr _1 (37) 
ei = (5 ±3) x 10 _15 s _1 (38) 
Xp = (-6.5 ± 0.8) x 10~ 13 (39) 

Only x p is clearly detected, while the two others are at best marginal detections. Note that proper 
motion can also lead to a change in the projected semi-major axis x p of the pulsar. However, if 
the observed x p were due to proper motion, Arzoumanian et al. (1996) find that the inner binary 
companion of the pulsar would then have a mass m c > 1.0 Mq, with sinii < 10° for a 1.35 
pulsar, which seems very unlikely. Hence we assume here that the observed x p is caused by the 
secular precession of the inner orbit induced by the presence of the second companion. 

To incorporate these measurements into our theoretical model, we have performed Monte- 
Carlo simulations with the unknown variables in the system, namely ii, 12, e 2 , 02, and 02- The 
angles 6*2 and 02 can be determined using i\, 12, w 2 , A 2 , and an additional undetermined angle 
a, which (along with i\ and 12) describes the relative orientation of the planes of the orbit 
of the inner binary and the orbit of the second companion. For a detailed description of the 
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geometry, see Appendix A. We assume a uniform probability distribution for cos(ii), cos (22), and 
a. Although there is no reason to expect the triple system to be in thermal equilibrium with the 
cluster (since the lifetime of the triple system in the core of M4 is only ~ 10 7 yr; see discussion in 
§3.3), we assume a thermal distribution for e2, i.e., a linear probability distribution Prob(e2)= 2e2 
(see, e.g., Heggie 1975) for lack of a better alternative. Our procedure for constructing random 
realizations of the system is as follows. We start by assuming a value of e2 and solving the 
nonlinear system numerically for mi-, a2, X2 and ll>2 as described in §3.1. Using this solution, we 
calculate r\ = (m2/mi)[(a p + a c )/ r i2] 3 - Then, for each trial, we generate random values for i±, 
12, and a, and calculate 62 and fa as described in Appendix A. We choose the number of trials 
for each assumed value of e2 so as to get a linear distribution for e2 over all the trials. We then 
calculate the secular perturbation rates using equations (33)-(36), and we check for consistency 
with the observed values of the perturbations. We use a simple rejection algorithm, accepting 
or rejecting a trial configuration based on a three-dimensional gaussian probability distribution 
which is the product of gaussian distributions for Cj\, e\ and x p centered around their mean values 
and with standard deviations equal to the error bars given above. 

In Figures 4 and 5, we show histograms of the number of successful trials for different values 
of 771,2 and the corresponding distance of the second companion from the inner binary ri2, for 
both f = f m (standard solutions) and /= 0.01 f m (required for obtaining stellar-mass solutions). 
For our standard solutions, we find that the mass distribution for the second companion peaks 
at 7772 = 0.01 ± 0.005 M , corresponding to a distance of ru = 38 ± 6 AU. Note that this is 
the distance of the second companion from the binary (not the semi-major axis of the second 
companion). Here the error bars represent 80% confidence intervals. The semi-major axis 02 
and period P2 are not well constrained and vary over three orders of magnitude. The period P2 
varies from about 10 2 to 10 4 years with the distribution centered around 400 years, and 02 varies 
from about 10 to 1000 AU with the distribution centered around 70 AU. In the second case, the 
distribution peaks at 7772 = 1.0 ± 0.5 M Q , corresponding to a distance of r\2 = 150 ± 35 AU. The 
period P2 varies from about 10 3 to 10 5 years with the distribution centered around 3000 years, 
and ci2 varies from about 10 2 to 10 4 AU with the distribution centered around 200 AU. We get 
peaks in the stellar mass range for 7772 only when —0.1 ^/ / f m £> 0.1. In all cases we find that 
the orientation of the second companion with respect to the inner binary (angles 62 and fa) is 
poorly constrained by the current data, as is the inclination 12 of the second companion. The 
inclination of the binary is better constrained to ~ 55° ± 15° (cf. Fig. 6). The eccentricity of the 
second companion, e2, is also poorly constrained. 

3.3. Model Predictions 

We can obtain a predicted value of the fifth pulse frequency derivative at the current 
epoch for each of our solutions by differentiating equation (9) for a specified value of e^. In 
Figure 7, we show the most probable values for from the Monte-Carlo simulations of §3.2. For 
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f = f m (standard solution), we find that the most probable value is = (0.15± 0.05) x 10 48 s 6 , 
whereas for /= 0.01 f m (which gives stellar mass solutions), /( 5 ) = -6.0 ± 0.4 x 10~ 51 s~ 6 . Thus 
even a crude measurement of /^5) should completely settle the question of the second companion's 
mass. 

We have also calculated the predicted evolution of the frequency derivatives / through / 
for the next 20 years. We show the results for the typical orbit of a Jupiter-to-brown-dwarf-size 
companion (m,2 = 0.01 M ) in Figure 8, and for a stellar-mass companion (777,2 = 0.5 M ) in 
Figure 9. In the first case, the orbit has e2 = 0.77, a period P2 = 1562 yr and 02 = 160 AU. We 
see that / changes sign in about 10 years, and that / decreases surprisingly fast, changing sign in 
about 1.5 years. For the stellar-mass case, only / changes sign in 10 years. The other frequency 
derivatives do not change significantly over 20 years. The orbit in this case has e2 = 0.49, a period 
P2 = 2034 yr and 02 = 161 AU. Thus a change in sign of / within the next couple of years would 
provide additional support for the existence of a planet or brown dwarf in this system. 

We also find that, in all cases, the values of /, /, and / at apastron are at least two, three, 
and five orders of magnitude smaller, respectively, than their present observed values. This means 
that the triple nature of the system would probably remain undetectable near apastron. It is 
therefore reasonable to find the second companion relatively close to periastron in our solutions 
(within ~ 15° for the case illustrated in Fig. 8, and ~ 40° for the case considered in Fig. 9). 

3.4. Discussion 

As mentioned previously, the method for determining the orbital parameters of a binary 
pulsar presented in §2, although quite general in its formulation, can only be applied successfully to 
systems containing fast millisecond pulsars in which the dynamically-induced frequency derivatives 
dominate the measurements, parts. In addition, it requires several successive higher-order 
frequency derivatives to be measured accurately. The PSR B1620— 26 triple system satisfies all 
of these conditions, and hence is ideally suited for analysis using our method. PSR B1620— 26 
has been observed for more than seven years, and its hierarchical triple structure is strongly 
supported by all current observations (Backer, Foster, & Sallmen 1993; Thorsett, Arzoumanian, 
&; Taylor 1993; Arzoumanian & Thorsett 1996). The error bar on / is likely to shrink rapidly as 
more timing data become available. If the actual value of / is close to the current best-fit value 
f m = —2.1 x 10~ 40 s -5 , then the second companion must have a mass 777.2 < 0.1 M as long as 
the system has an inclination 12 > 7°, with the most probable mass (given by our Monte-Carlo 
simulations) being 0.01 ± 0.005 M (the error-bar indicates an 80% confidence interval). If / is 
within la of f m , then the same result holds to within a factor of two. A rather low inclination 
angle {12 < 10°) or | / / f m \ < 0.1 (i.e., / more than 3<r away from f m ) would be required if 
the second companion is a main-sequence star with 777-2 > 0.1 M . 

Instead our results clearly suggest that the second companion is a ~ 0.01 M brown dwarf 
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or giant planet. This is surprising since low-mass objects are not expected to be found in the 
cores of globular clusters. The reason is that low-mass objects have higher velocities in energy 
equipartition and are preferentially ejected from globular clusters as they evaporate in the tidal 
field of the Galaxy. Hence we do not expect to find very low-mass stars or brown dwarfs in 
globular clusters, especially not near the core. Recent HST observations of globular clusters 
(e.g., Paresce, De Marchi, & Romaniello 1995) also support this view by finding that stellar mass 
functions in clusters flatten or even drop for masses below ~ 0.1 M . In addition, if the second 
companion of PSR B1620— 26 is indeed of low mass, then the unusually high eccentricity of the 
inner binary pulsar cannot be explained by secular perturbations due to the second companion, 
since that would require a stellar-mass second companion (Rasio 1994). It would also preclude any 
possibility of an optical identification of the triple system. Bailyn et al. (1994) have searched deep 
optical images of M4 for an optical counterpart of the pulsar. They have identified a candidate 
which, if interpreted as a single object, could be a 0.45 M main-sequence star within 0.3" of the 
nominal pulsar position. However, it is possible that this object is in fact a blend of fainter stars 
not associated with the pulsar, or simply a chance superposition. Future observations of the region 
with HST, as well as improved ground-based astrometry, should help resolve the issue. 

Low-mass stars and brown dwarfs could exist in dense globular cluster cores as binary 
companions to more massive stars. Dynamical interactions could then lead to an exchange, leaving 
the low-mass object in orbit around a neutron star. Indeed, Sigurdsson (1992) had discussed the 
possibility of finding planetary companions to pulsars in globular clusters even before the triple 
nature of PSR B1620— 26 was established. A possible formation scenario for the triple system 
starts with an interaction between a neutron-star-white-dwarf binary and a main-sequence star 
with a large Jupiter-type planet or brown-dwarf companion (cf. Sigurdsson 1995, 1993). As a 
result of this interaction the white dwarf is ejected while the main-sequence star and its planet 
or brown-dwarf companion remain in orbit around the neutron star. The main-sequence star, 
as it evolves and later expands as a red giant, would then transfer mass onto the neutron star, 
thus spinning it up and forming the millisecond pulsar in the triple configuration we see today. 
However, tidal dissipation during the mass transfer phase would effectively circularize the orbit 
of the binary, leaving a residual eccentricity e± Ss 10 -4 (Phinney 1992). Therefore this formation 
scenario leaves the much higher observed eccentricity (ei ~ 0.03) of the inner binary unexplained. 
It has been suggested that the eccentricity of the inner binary may have been induced during a 
dynamical interaction with another cluster star. The probability of disrupting the triple during 
such an interaction is only ~ 0.5 (Sigurdsson 1995). Based on the results of Rasio & Heggie (1995), 
however, we find that the observed eccentricity would require an encounter with a distance of 
closest approach of ~ 2.5 AU, considerably smaller than the size of the outer orbit, and occuring on 
average once in ~ 4 x 10 8 yr. For comparison, the lifetime of the triple system in the cluster is only 
about r ~ 10 8 yr/?4 1 (T5 (02/10 AU)^ 1 ~ 2 x 10 7 yr, where p = 10 4 ^4 M pc~ 3 is the density near 
the center of M4, a = ba^kms' 1 is the velocity dispersion, and ci2 is the size of the outer orbit 
(Rasio 1994). For one interaction that could have produced the eccentricity of the inner binary, 
we therefore expect ~ 20 interactions that could have disrupted the triple, each with probability 
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~ 0.5, leaving the probability of survival at ~ 10~ 6 . An additional problem is that the age of the 
millisecond pulsar in this scenario must be comparable to the age of the triple (~ 10 7 yr), which 
requires the millisecond pulsar to be extremely young. This problem could be avoided if the triple 
was instead formed during an interaction involving a pre-existing binary millisecond pulsar and 
another primordial binary (containing the present second companion and another star that was 
ejected during the interaction; see Rasio, McMillan, & Hut 1995). The current eccentricity of 
the (inner) binary pulsar could then have been induced during the same interaction that formed 
the triple, although this would require some fine tuning. More significantly, one would expect 
the more massive member of the other binary rather than the low-mass object (Jupiter or brown 
dwarf) to be preferentially retained in the triple while the other gets ejected. 

Naturally, if the low-mass object (Jupiter or brown dwarf) was attached to a much more 
massive star, it is easier to understand how it was retained by the cluster and why it is now 
found close to the cluster core. In particular, in the first formation scenario discussed above, 
the main-sequence star must have been fairly massive (m ~ 1M Q ) to have evolved into a red 
giant after the triple was formed. Two-body relaxation in the cluster will tend to bring this 
main-sequence star (with its attached low-mass companion) down to the cluster core since it is 
more massive than the average object in the cluster. Confirmation of the existence of a ~ 0.01 M Q 
object in PSR B1620— 26 would therefore provide further indication that many stars, even in 
globular clusters, could have very low-mass companions or planets. This is especially important 
in the light of recent discoveries of several ~ 10 -3 — 10 -2 M Q objects around nearby stars (e.g., 
Mayor & Queloz 1995; Butler & Marcy 1996; Marcy k Butler 1996). 

4. Application to the PSR B1257+12 Planetary System 

We now turn to the application of our method to the planetary system around the millisecond 
pulsar PSR B1257+12. This system contains three confirmed Earth-mass planets in quasi-circular 
orbits (Wolszczan & Frail 1992; Wolszczan 1994). The planets have masses of 0.015/ siniiM^, 
3.4/ sini2 M ffi , and 2.8/ sin 23 M ffi , where ii, 22 and i^ are the inclinations of the orbits with respect 
to the line of sight, and are at distances of 0.19 AU, 0.36 AU, and 0.47 AU, respectively, from 
the pulsar. In addition, the unusually large second and third frequency derivatives of the pulsar 
suggest the existence of a fourth, more distant and massive planet in the system (Wolszczan 1996). 

4.1. Analysis of the Frequency Derivative Data 

The residual pulse frequency derivatives for PSR B1257+12 (after subtraction of a model 
for the inner three planets) are / = -8.6 x 10" 16 s~ 2 , / = (-1.25 ± 0.05) x 10~ 25 s" 3 , 
and f= (1.1 ± 0.3) x 10~ 33 s~ 4 (Wolszczan 1996), while the frequency / = 160.8s- 1 . The 
value of / has been corrected for the apparent acceleration due to the pulsar's transverse 
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velocity (the so-called Shlovskii effect; cf. Camilo et al. 1994). The error bars on / and /, 
for the purposes of this discussion, are negligible. Note that the measurement of / is only 
preliminary, but we assume here that the value quoted above (from Wolszczan 1996) is correct. 
Comparison with PSR B1855+09, which has a very similar pulse frequency, / = 186 s _1 , and first 
frequency derivative / = —6.2 x 10~ 16 s~ 2 (Kaspi et al. 1994) indicates that the observed / for 
PSR B1257+12 could well be entirely (or in large part) intrinsic rather than acceleration-induced. 
The timing age for the pulsar, r = —f/2f ~3x 10 9 yr, is entirely consistent with that expected 
for a millisecond pulsar. Therefore we will treat / acc essentially as a free parameter in our analysis. 
The observed /, on the other hand, is two orders of magnitude larger than for PSR B1855+09, 
which has / < 2.0 x 10 -27 s~ 2 . Thus the observed / is almost certainly due to the presence of 
another planet rather than intrinsic timing noise in the pulsar. 

With three frequency derivatives measured, we can use the method of §2.3 to model the 
system. Given the near-circular orbits of the three inner planets, it is natural to assume that 
the orbit of the fourth planet also has a low eccentricity. In addition, it is easy to show that 
dynamical interactions with passing stars in the Galaxy are not likely to produce any significant 
perturbations of the system (which could otherwise increase the eccentricity of an outer planet's 
orbit; cf. Heggie & Rasio 1996). 

Since the value of / acc is uncertain, we explore a wide range, 0.01 < / acc //obs < 1- Note that, 
for a circular orbit, / acc and / must have opposite signs (cf. eqs. [19] and [21]). Hence / acc cannot 
be positive. For each value of / acc , we calculate the mass and semi-major axis of the fourth planet 
using equations (25), (28) and (30). We illustrate the results in Figure 10. We find that the mass 
of the fourth planet varies significantly, from ~ 0.08 M® (for / acc = 0.01/ o b s ) to ~ 100 M© (for 
/acc = /obs)- The simplest interpretation of the present best-fit values of the frequency derivatives, 
assuming / acc = / b s > implies a mass of about 100/ sin 24 M® (i.e., comparable to Saturn's mass) 
for the fourth planet, at a distance of about 38 AU (i.e., comparable to Pluto's distance from the 
Sun), and with a period of about 170 yr in a circular, coplanar orbit (Wolszczan 1996). However, 
if /acc 7^ /obs, then the fourth planet can have a wide range of masses. In particular, it can have 
a mass comparable to that of Mars (at a distance of 9 AU), Uranus (at a distance of 25 AU) or 
Neptune (at a distance of 26 AU), for / acc = 0.015, 0.30, or 0.34 / Q b s , respectively. 



4.2. Discussion 

In this system the perturbations of the inner planets produced by the fourth planet are 
probably far too small to be detected. This is in contrast to the mutual perturbations of the inner 
planets themselves, which are important and have been detected (Rasio et al. 1992; Wolszczan 
1994). Using equations (33)-(36), we predict e ~ 10 - s -1 and Co ~ 10~ 7 degyr _1 for the orbit of 
the third planet, assuming that all orbits are coplanar and that the mass of the fourth planet is 
100 M®. The perturbations for the two innermost planets are even smaller. Hence the existence of 
the fourth planet is likely to be confirmed only through further measurements of pulse frequency 
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derivatives. 

It has been pointed out that the masses and radii of the three inner planets in PSR B1257+12 
are in the same ratios as the masses and radii of the corresponding first three planets in the 
Solar System (Mazeh & Goldman 1995). This might perhaps be indicative of a global underlying 
formation mechanism for the two systems. 

Although the fourth planet could have the same mass (normalized to the mass of the third 
planet) as that of Mars, Uranus, or Neptune (normalized to the mass of the Earth), the ratio of 
radii in each case would be much larger than the corresponding ratio for the solar system (cf. 
Fig. 10). Thus this system does not seem to maintain its regularity with the Solar System, since 
the mass and radius ratios of the fourth planet would not simultaneously match those of any 
planet in the Solar System. This is true for the entire range of values of / acc considered above. 
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A. Geometry of the Triple Configuration 

The orbit of the inner binary and the orbit of the second companion in general do not lie in 
the same plane. The inclinations of the two planes with respect to the line of sight are given by 
i\ and i^. To specify the plane of an orbit completely, one needs the inclination angle together 
with another azimuthal angle a (which lies between and 2ir). Since the reference axis for a is 
arbitrary, we can take it to lie in the plane of one of the orbits, so that a is the difference between 
the azimuthal angles of the two planes. In random Monte-Carlo trials, a is then taken to be 
uniformly distributed between and 2ir. 

In order to determine 02 and <p2 using the other angles, we need to change coordinates between 
two reference frames. The first frame has its origin at the center of mass of the inner binary, 
with the x-axis in the plane of the orbit of the second companion, the y-axis passing through the 
pericenter of the orbit, and the z-axis perpendicular to the plane of the orbit, so that the motion of 
second companion is counterclockwise around the z-axis. We shall refer to the coordinates of the 
second companion in this frame as (x' c , y' c , z' c ). Then, x' c = — ri2sinA2, y' c = ri2CosA2, and z' c = 0. 

The second frame similarly has its origin at the center of mass of the inner binary, with the 
x-axis in the plane of the orbit of the inner binary, the y-axis passing through the pericenter of 
the orbit, and the z-axis perpendicular to the plane of the orbit, so that the motion of pulsar is 
counterclockwise around the z-axis. We wish to find the coordinates (x", y", z") of the second 
companion in this frame. We can then calculate 62 and <j>2 using the formulae 62 = cos -1 ^"/?"^) 
and 4>2 = tan _1 (— x"/y"), keeping in mind that for y" c < 0, we must add 180° to <p2 in order to get 
the correct quadrant. 

In order to get (x", y", z") from (x' c , y' c , z' c ), we shall rotate the first frame to the second 
frame using the standard Euler angles formalism (see, e.g., Goldstein 1980). In this formalism, 
any arbitrary rotation of an object is represented as a sequence of three consecutive rotations- first 
about the z-axis by an angle (p, then about the new x-axis by an angle 6, and finally again about 
the new z-axis by an angle ip. 

In order to use this formalism, we will use an intermediate frame of reference that is fixed in 
space, with its origin at the center of mass of the inner binary, the y-axis along the line of sight, 
the x-axis in the plane of the orbit of the second companion, and the z-axis such that the motion 
of the second companion is counterclockwise about the z-axis. The first frame described above 
can be obtained from this fixed frame by rotating it through the Euler angles 0, %2 and u>2 — 90°, 
respectively. Similarly, the second frame can be obtained from the fixed frame by rotating it 
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through the angle a about the y-axis, and then rotating it through the Euler angles 0, i\ and 
Ui — 90°, respectively. The sequence of Euler angle rotations is represented as a matrix A(cp, 9, ip) 
We get the coordinates (x", y'^ z") from (x' c , y' c , z' c ) by multiplying first by the inverse matrix 
A _1 (0, i2,uJ2 — 90°), then by multiplying by the matrix for rotation about the y-axis B(a), and 
then multiplying by the matrix A(0, h,u>i — 90°). The matrices are given by, 



A(^0,y>) 



/ cos ip cos 4> — cos 9 sin <p sin ip cos ip sin cp + cos 9 cos <p sin ip sin ip sin 9 \ 

— sin ip cos (f) — cos 9 sin cos tp — sin ^ sin + cos 9 cos cos tp cos ^ sin 9 
\ sin 6* sin — sin 9 cos cos 9 J 



and 



B(a) 



/ cos a — sin a \ 

1 
\ sina cos a / 
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Fig. 1. — Allowed values of the semi-major axis 02, mass 777-2, eccentricity e2, longitude at epoch 
A2, and longitude of pericenter u>2 for the second companion of PSR B1620— 26, using the latest 
available values for all pulse frequency derivatives. This is our "standard solution," using the 
present best-fit value f m = —2.1 x 10 -40 s -5 and assuming that all measured frequency derivatives 
are dynamically induced. Acceptable solutions all have 2.4 x 10 -4 M Q < 1112 sini2 < 1.2 x 10 -2 M Q . 
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Fig. 2. — The orbital period, mass and eccentricity of the second companion of PSR B1620— 26, 
for different values of the acceleration-induced first frequency derivative. The various curves arc 
for / acc = / obs (solid line), / acc = 0.1 / o6s (long-dashed line), and f acc = 0Mf obs (short-dashed 
line). We assume that the inclination angle %2 = 90°. We see that the mass range does not change 
significantly when varying / acc , and stellar masses are always excluded. 
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Fig. 3. — Similar to Fig. 2, but the dependence of the solutions on / is illustrated for an expanded 
4(7 error bar around the best-fit value f m = —2.1 x 10~ 40 s~ 5 . The solutions are shown for f = f m 
(solid line), f = f m —la (long-dashed line), f = f m +lc (short-dashed line), /= 0.1 f rn (dot- 
dashed line), and /= 0.01 f m (dotted line), where a = 0.6 x 10~ 40 s~ 5 . Stellar-mass solutions are 
obtained when —0.1 / f m £> 0.1 
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Fig. 4. — Histogram of the number of successful trials (N) for different values of rri2 and the 
corresponding distance of the second companion from the inner binary ri2, for the case f = f m 
in our Monte-Carlo simulations. We find that the most probable value for the second companion 
mass is iri2 = 0.01 ± 0.005 M Q , corresponding to a distance of ru = 38 ± 6AU (80% confidence 
intervals) . 
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Fig. 5. — Same as Fig. 4, but for the case / = 0.01 f m . We see that the most probable value for the 
second companion mass is now ttt-2 = 1.0±0.5 M , corresponding to a distance of m = 150±35 AU. 
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Fig. 6. — Histograms of the number of successful trials (N) for various parameters in the Monte- 
Carlo simulations with f = f m . All the angles are poorly constrained except the inclination of the 
inner binary, which is slightly better constrained to i\ ~ (55 ± 15)°. 
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Fig. 7. — The most probable value of the fifth frequency derivative given by the Monte-Carlo 
simulations. For f = f m (lower frame), /( 5 ) ss 0.15(5) x 10~ 48 s~ 6 and for /= 0.01 f m (upper 
frame), /( 5 ) « -0.0060(4) x 10~ 48 s- 6 . 
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Fig. 8. — The predicted variation of the frequency derivatives over the next 20 years, for a low- 
mass (777.2 = 0.01 Mq) second companion (standard solution). The units for the four derivatives are 
10~ 15 s~ 2 for /, 10 -23 s~ 3 for /, 10~ 33 s -4 for /, and 10~ 40 s" 5 for /. We see that / changes sign 
in about 10 years, and / changes sign in about 1.5 years. 
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Fig. 9. — The predicted variation of the frequency derivatives over the next 20 years, for a stellar- 
mass (rrt2 = 0.5 Mq) second companion (assuming /= 0.01 f m ). The units for the four derivatives 
are 10~ 15 s~ 2 for /, 10~ 23 s -3 for /, 10~ 33 s~ 4 for /, and 10~ 42 s -5 for /. We see that / changes 
sign in about 10 years, but the other derivatives do not change much in 20 years. 
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Fig. 10. — The mass and semi-major axis of the possible outer planet in the PSR B1257+12 
planetary system for / acc in the range 0.01-1.0 / bs- The present best-fit values of the frequency 
derivatives with / acc = / Q b s imply the presence of a planet with mass ~ 100/ sin 14 M®, at a distance 
of « 38 AU. The marked points on the curve indicate the values of face/ fobs- The points labelled 
M, U and N indicate configurations with the same mass and radius ratios (in this system) as those 
of Mars, Uranus and Neptune (in the Solar System), respectively. 



